pacman::p_load(SpatialAcc, sf, tidyverse, tmap, ggstatsplot)
In-Class Exercise 9
In this exercise, we will import geospatial data, create buffers for accessibility analysis, use spatial joins to assess facility proximity, visualize spatial data on maps, and calculate Hansen’s Accessibility metrics.
1 Exercise Reference
2 Overview
In this exercise, we will import geospatial data, create buffers for accessibility analysis, use spatial joins to assess facility proximity, visualize spatial data on maps, and calculate Hansen’s Accessibility metrics.
3 Import the R Packages
The following R packages will be used in this exercise:
| Package | Purpose | Use Case in Exercise |
|---|---|---|
| tidyverse | Data manipulation and visualization in R. | Data wrangling, cleaning, and visualization. |
| sf | Handling and visualizing geospatial data. | Importing and managing geospatial data, creating buffers. |
| SpatialAcc | Measuring spatial accessibility metrics. | Computing Hansen’s Accessibility. |
| tmap | Thematic map visualization. | Creating maps for spatial data visualization. |
| ggstatsplot | Enhanced statistical data visualization. | Visualizing accessibility differences across regions. |
To install and load these packages, use the following code:
4 Data Wrangling
We will download and use the CHAS Clinics and Eldercare Services data sets from data.gov.sg portal.
To read the ELDERCARE shapefile data:
eldercare <- st_read(dsn = "data/raw_data",
layer = "ELDERCARE") %>%
st_transform(crs = 3414)
Reading layer `ELDERCARE' from data source
`/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/raw_data'
using driver `ESRI Shapefile'
Simple feature collection with 133 features and 18 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21
To read the CHAS clinics kml file:
CHAS <- st_read("data/raw_data/CHASClinics.kml") %>%
st_transform(crs = 3414)
Reading layer `MOH_CHAS_CLINICS' from data source
`/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/raw_data/CHASClinics.kml'
using driver `KML'
Simple feature collection with 1193 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84
We will also use the datasets from Hands-on Exercise 9
mpsz <- st_read(dsn = "data/geospatial",
layer = "MP14_SUBZONE_NO_SEA_PL") %>%
st_transform(crs = 3414)
Reading layer `MP14_SUBZONE_NO_SEA_PL' from data source
`/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
hexagons <- st_read(dsn = "data/geospatial",
layer = "hexagons") %>%
st_transform(crs = 3414)
Reading layer `hexagons' from data source
`/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 3125 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 21506.33 xmax: 50010.26 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
eldercare <- st_read(dsn = "data/geospatial",
layer = "ELDERCARE") %>%
st_transform(csr = 3414)
Reading layer `ELDERCARE' from data source
`/Users/walter/code/isss626/isss626-gaa/In-class_Ex/In-class_Ex09/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 120 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
Projected CRS: SVY21 / Singapore TM
ODMatrix <- read_csv("data/aspatial/OD_Matrix.csv",
skip = 0)
4.1 Data Cleaning and Updating Attributes
4.1.1 Supply:
Set to 100 first to calibrate the model.
eldercare <- eldercare %>%
select(fid, ADDRESSPOS) %>%
mutate(capacity = 100)
4.1.2 Demand:
hexagons <- hexagons %>%
select(fid) %>%
mutate(demand = 100)
4.1.3 OD Matrix
Convert the distance matrix to km units:
distmat_km <- as.matrix(distmat/1000)
5 Count Number of CHAS Clinics within a Distance Around Each Eldercare Centre
To count the number of points within a distance, we will do the following steps:
- create a buffer of 1 km around each eldercare centre
- visualize the data
- count points
5.1 Create Buffer
Note the singapore crs is in metres. To create a 1km buffer,
dist should be 1000.
buffer_1km <- st_buffer(eldercare,
dist = 1000)
5.2 Visualize Plots
tmap_mode("view")
tm_shape(buffer_1km) +
tm_polygons() +
tm_shape(CHAS) +
tm_dots()
5.3 Count Points
To count number of CHAS Clinics within 1 KM around each Eldercare Centre
buffer_1km$pts_count <- lengths(
st_intersects(buffer_1km, CHAS))
buffer_1km$pts_count
[1] 17 31 26 13 26 28 13 24 11 14 17 32 6 20 21 32 15 19 23 19 11 17 20 22 41
[26] 17 18 39 24 25 15 16 28 28 14 30 23 6 21 15 9 24 24 25 30 8 32 22 17 28
[51] 33 16 16 16 8 23 16 12 17 19 8 21 14 20 12 22 11 24 18 17 13 16 13 18 19
[76] 12 19 23 23 21 15 10 15 23 9 22 26 28 22 16 14 31 25 30 9 14 14 14 10 22
[101] 34 19 16 11 9 17 24 12 31 28 17 16 26 31 24 32 22 29 18 30
6 Computing Hansen’s Accessibility
To compute Hansen’s Accessibility:
acc_Hansen <- data.frame(ac(hexagons$demand,
eldercare$capacity,
distmat_km,
#d0 = 50,
power = 2,
family = "Hansen"))
Then we tidy the output:
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hexagon_Hansen <- bind_cols(hexagons, acc_Hansen)
6.1 Visualising Hansen’s Accessibility
mapex <- st_bbox(hexagons)
tmap_mode("plot")
tm_shape(hexagon_Hansen,
bbox = mapex) +
tm_fill(col = "accHansen",
n = 10,
style = "quantile",
border.col = "black",
border.lwd = 1) +
tm_shape(eldercare) +
tm_symbols(size = 0.1) +
tm_layout(main.title = "Accessibility to eldercare: Hansen method",
main.title.position = "center",
main.title.size = 2,
legend.outside = FALSE,
legend.height = 0.45,
legend.width = 3.0,
legend.format = list(digits = 6),
legend.position = c("right", "top"),
frame = TRUE) +
tm_compass(type="8star", size = 2) +
tm_scale_bar(width = 0.15) +
tm_grid(lwd = 0.1, alpha = 0.5)
6.2 Visualize Statistical Graphic
hexagon_Hansen <- st_join(hexagon_Hansen, mpsz,
join = st_intersects)
ggbetweenstats(
data = hexagon_Hansen,
x = REGION_N,
y = accHansen,
type = "p")